Demographics and reproductive biology of H. schistosus from bycatch on the west coast of India
Analysis for Manuscript
Introduction
This the code for analysis and supplementary figures for the manuscript. The manuscript is currently under review.
#Required libraries
require(tidyverse)
require(lubridate)
require(viridis)
# Importing data
snakes = read.csv("./Data/Sea-snakes_fish-dep_mastersheet_250420.csv")
embryos = read.csv("./Data/Sea snake_embryo data_mastersheet.csv")
snakes <- snakes%>%
mutate(Year = as.factor(Year))%>%
filter(Species == "Hydrophis curtus" | Species == "Hydrophis schistosus")%>%
droplevels()
HS <- snakes%>%
filter(Species == "Hydrophis schistosus")%>% # filtering out other snakes from bycatch data
mutate(Year = as.factor(Year))
# function to calculate mcfadden's r2
mcfadden <- function(x){
r2 <- 1 - (x$deviance/x$null.deviance)
}
# sampling years
years <- data.frame(Year = factor(c(2016:2019)))Sample size
The number of snakes (N) sampled from bycatch of trawlers, shore seines and gillnetss.
snakes%>%
group_by(Species)%>%
count(name = "N")| Species | N |
|---|---|
| Hydrophis curtus | 282 |
| Hydrophis schistosus | 967 |
The number of snakes sampled in each year:
n.yr <- snakes%>%
group_by(Species, Year)%>%
count(name = "N")
n.yr%>%
spread(Year, N)| Species | 2016 | 2017 | 2018 | 2019 |
|---|---|---|---|---|
| Hydrophis curtus | 75 | 38 | 123 | 46 |
| Hydrophis schistosus | 178 | 215 | 521 | 53 |
The number of trips (n) and fishing effort (haul-hours) sampled for sea snakes in the current study.
snakes%>%
filter(Gear.Type != "")%>%
select(Date, Gear.Type, No..of.Hauls, Average.Haul.Duration..Hours., Tow.duration.hours)%>%
distinct()%>%
# Calculating tow duration
mutate(Tow.duration.hours = case_when(!is.na(Tow.duration.hours) ~ Tow.duration.hours,
is.na(Tow.duration.hours) ~ No..of.Hauls*Average.Haul.Duration..Hours.))%>%
group_by(Gear.Type)%>%
summarise(n = n(), # counting number of trips sampled
haul.hours = sum(Tow.duration.hours, na.rm = T)) # total fishing effort| Gear.Type | n | haul.hours |
|---|---|---|
| GillNet | 340 | 535.9633 |
| Rampan | 46 | 190.4500 |
| Trawler | 76 | 285.1600 |
Demographics of H. schistosus across time
Age structure
# calculating the mean SVL across years
age.yr <- snakes%>%
group_by(Species, Year)%>%
summarise(N = n(),
mean = mean(Snout.to.Vent..cm., na.rm = T))
# SVL cut off for age classes
maturtity <- snakes%>%
group_by(Species)%>%
count()%>%
mutate(juv = 35,
adult = ifelse(Species == "Hydrophis curtus", 54, 65))Distribution fo SVL across years
# plotting distribution of SVL across years
snakes%>%
filter(Snout.to.Vent..cm. > 20)%>%
ggplot(aes(Year, Snout.to.Vent..cm.))+
geom_violin(fill = "grey")+
geom_boxplot(width = 0.1)+
geom_hline(data = maturtity, aes(yintercept = adult), linetype = "dashed")+
geom_hline(data = maturtity, aes(yintercept = juv), linetype = "dotted")+
#stat_summary(fun.data = "mean_sdl", geom = "pointrange", size = 1)+
geom_label(data = age.yr, aes(Year, 10, label = N))+
facet_wrap(~Species, ncol = 1, scales = "free_y")+
labs(y = "Snout to vent length (cm)")+
theme(strip.text = element_text(face = "italic"))ggsave(last_plot(), filename = "./Figures/figure1.tiff", height = 6, width = 8)We did not find any H. curtus neonates in our sampling. Majority of H. curtus were juveniles. H. schistosus populations consisted of mostly adults.
Linear model for change in SVL of each species across years:
library(car)
snakes%>%
group_by(Species)%>%
select(Year, Snout.to.Vent..cm.)%>%
nest()%>%
mutate(mod = map(data, ~lm(Snout.to.Vent..cm. ~ Year, data = .)),
sumr = map(mod, broom::tidy),
stat = map(mod, glance))%>%
select(sumr, stat)%>%
unnest()%>%
select(c(Species:r.squared))| Species | term | estimate | std.error | statistic | p.value | r.squared |
|---|---|---|---|---|---|---|
| Hydrophis schistosus | (Intercept) | 85.598303 | 1.850695 | 46.2519793 | 0.0000000 | 0.0689121 |
| Hydrophis schistosus | Year2017 | -13.996389 | 2.122893 | -6.5930743 | 0.0000000 | 0.0689121 |
| Hydrophis schistosus | Year2018 | -10.158177 | 1.974836 | -5.1438083 | 0.0000003 | 0.0689121 |
| Hydrophis schistosus | Year2019 | -1.732303 | 2.818894 | -0.6145327 | 0.5390387 | 0.0689121 |
| Hydrophis curtus | (Intercept) | 65.043404 | 1.511818 | 43.0233133 | 0.0000000 | 0.1470604 |
| Hydrophis curtus | Year2017 | -5.212154 | 2.405287 | -2.1669571 | 0.0312836 | 0.1470604 |
| Hydrophis curtus | Year2018 | -8.833123 | 1.825449 | -4.8388759 | 0.0000024 | 0.1470604 |
| Hydrophis curtus | Year2019 | -12.998166 | 2.225335 | -5.8409920 | 0.0000000 | 0.1470604 |
ANOVA to test for effect of years:
snakes%>%
group_by(Species)%>%
select(Year, Snout.to.Vent..cm.)%>%
nest()%>%
mutate(mod = map(data, ~lm(Snout.to.Vent..cm. ~ Year, data = .)),
stat = map(mod, car::Anova))%>%
select(stat)%>%
unnest()| Species | Sum Sq | Df | F value | Pr(>F) |
|---|---|---|---|---|
| Hydrophis schistosus | 13334.508 | 3 | 19.66265 | 0e+00 |
| Hydrophis schistosus | 180165.620 | 797 | NA | NA |
| Hydrophis curtus | 4363.959 | 3 | 12.98867 | 1e-07 |
| Hydrophis curtus | 25310.652 | 226 | NA | NA |
Plotting fitted values:
snakes%>%
group_by(Species)%>%
select(Year, Snout.to.Vent..cm.)%>%
nest()%>%
mutate(mod = map(data, ~lm(Snout.to.Vent..cm. ~ Year, data = .)),
fit = map(mod, ~predict(., years, se.fit = T)),
fit = map(fit, ~as.data.frame(.)))%>%
select(fit)%>%
unnest()%>%
mutate(Year = factor(c(2016:2019)))%>%
ggplot(aes(Year, fit))+
geom_pointrange(aes(ymin = fit - se.fit, ymax = fit + se.fit))+
labs(y = "SVL (cm)")+
facet_wrap(~Species)+
theme(strip.text = element_text(face = "italic"))There was a significant change in SVL distribution of species across years.
Proportion of developmental classes across years
pop.yr <- snakes%>%
filter(!is.na(Class))%>% # fix 2016 svl
group_by(Species, Year)%>%
# Counting number of individuals in each developmental class
count(Class)%>%
complete(Class, fill = list(n = 0))%>%
# Total number sampled in each year
mutate(N = sum(n))%>%
group_by(Year, Species, Class)%>%
# Proportion of each developmental class
mutate(prop.age = n/N)
# Summary stats
pop.yr%>%
group_by(Species, Class)%>%
skimr::skim(prop.age)%>%
skimr::yank("numeric")%>%
select(c(Species, Class, mean, sd))Variable type: numeric
| Species | Class | mean | sd |
|---|---|---|---|
| Hydrophis curtus | Adult | 0.43 | 0.31 |
| Hydrophis curtus | Juvenile | 0.57 | 0.31 |
| Hydrophis curtus | Neonate | 0.02 | NA |
| Hydrophis schistosus | Adult | 0.88 | 0.09 |
| Hydrophis schistosus | Juvenile | 0.14 | 0.07 |
| Hydrophis schistosus | Neonate | 0.03 | 0.02 |
Testing change in proportion of developmental classes across years
Fitting a multinomial logistic model to determine change in age classes across years in each species.
require(nnet)
require(broom)
dev.yr <- snakes%>%
select(Species, Year, Class)%>%
filter(!is.na(Class))%>% # removing empty class
group_by(Species)%>%
nest()%>%
mutate(mod = map(data, ~multinom(Class ~ Year, data = .)),
null = map(data, ~multinom(Class ~ 1, data = .)),
smry = map(mod, ~broom::tidy(., exonentiate = T)),
# Calculating mc fadden's r2
r2 = map2_dbl(mod, null, ~ 1 - .x$value/.y$value))dev.yr%>%
select(smry, r2)%>%
unnest()| Species | y.level | term | estimate | std.error | statistic | p.value | r2 |
|---|---|---|---|---|---|---|---|
| Hydrophis schistosus | Juvenile | (Intercept) | -2.4450606 | 0.7380571 | -3.3128339 | 0.0009236 | 0.0455695 |
| Hydrophis schistosus | Juvenile | Year2017 | 1.1234010 | 0.7573074 | 1.4834147 | 0.1379643 | 0.0455695 |
| Hydrophis schistosus | Juvenile | Year2018 | 0.5064641 | 0.7510767 | 0.6743174 | 0.5001096 | 0.0455695 |
| Hydrophis schistosus | Juvenile | Year2019 | -13.3315259 | 377.0158759 | -0.0353606 | 0.9717922 | 0.0455695 |
| Hydrophis schistosus | Neonate | (Intercept) | -3.1352560 | 1.0212895 | -3.0698995 | 0.0021413 | 0.0455695 |
| Hydrophis schistosus | Neonate | Year2017 | -14.5022162 | 526.2494814 | -0.0275577 | 0.9780149 | 0.0455695 |
| Hydrophis schistosus | Neonate | Year2018 | -0.9347455 | 1.0900891 | -0.8574946 | 0.3911716 | 0.0455695 |
| Hydrophis schistosus | Neonate | Year2019 | -8.7344420 | 53.4649902 | -0.1633675 | 0.8702291 | 0.0455695 |
| Hydrophis curtus | Juvenile | (Intercept) | -1.9386218 | 1.0661369 | -1.8183610 | 0.0690090 | 0.0724193 |
| Hydrophis curtus | Juvenile | Year2017 | 2.7268690 | 1.1322946 | 2.4082681 | 0.0160284 | 0.0724193 |
| Hydrophis curtus | Juvenile | Year2018 | 2.5058458 | 1.0853066 | 2.3088829 | 0.0209501 | 0.0724193 |
| Hydrophis curtus | Juvenile | Year2019 | 3.5480223 | 1.1437099 | 3.1022049 | 0.0019208 | 0.0724193 |
| Hydrophis curtus | Neonate | (Intercept) | -11.1507510 | 99.7723289 | -0.1117620 | 0.9110121 | 0.0724193 |
| Hydrophis curtus | Neonate | Year2017 | -2.2749071 | 278.6764080 | -0.0081633 | 0.9934867 | 0.0724193 |
| Hydrophis curtus | Neonate | Year2018 | 8.2062502 | 99.7749668 | 0.0822476 | 0.9344498 | 0.0724193 |
| Hydrophis curtus | Neonate | Year2019 | -3.4778373 | 576.2501607 | -0.0060353 | 0.9951846 | 0.0724193 |
Plotting fitted estimates of proportions:
dev.yr%>%
mutate(fit = map(mod, ~data.frame(predict(., years, "probs"))))%>%
select(fit)%>%
unnest()%>%
mutate(Year = factor(c(2016:2019)))%>% # adding years variable
left_join(n.yr, by = c("Species", "Year"))%>% # getting sample size in each year
gather(c(Adult, Juvenile, Neonate), key = "Class", value = "p.hat")%>%
mutate(se = sqrt(p.hat*(1-p.hat)/N))%>% # caluclating standard error
# Plotting
ggplot(aes(Year, p.hat, color = Species))+
geom_pointrange(aes(ymin = p.hat-se, ymax = p.hat+se))+
labs(y = "Estimated proportion")+
facet_wrap(~Class)+
theme(legend.text = element_text(face = "italic"))Age structure of H. schistosus does not change significantly over a four year period from 2016 to 2018.
The proportion of H. curtus juveniles increases significantly between 2016 and 2018
Change in SVL distribution over months
month.svl <- HS%>%
filter(Year == "2018", # data not sufficient for other years
Snout.to.Vent..cm. > 20)%>% # removing erroneous data
mutate(Month = factor(Month, levels = month.name))%>%
complete(Month)%>%
group_by(Month)%>%
# Calculating mean SVL in each month
summarise(mean.SVL = mean(Snout.to.Vent..cm., na.rm = T))
# Month of observed births
births <- data.frame(Species = "Hydrophis schistosus", Month = "April")# plotting distribution of SVL across months
HS%>%
filter(Year == "2018", # data not sufficient for other years
Snout.to.Vent..cm. > 20)%>% # removing erroneous data
mutate(Month = factor(Month, levels = month.name))%>%
complete(nesting(Species), Month)%>%
droplevels()%>%
ggplot(aes(Month, Snout.to.Vent..cm.))+
geom_violin(fill = "light grey")+
#geom_point(data = month.svl, aes(x = Month, y = mean.SVL), size = 3)+
geom_boxplot(width = 0.1)+
geom_segment(data = births, aes(x = Month, xend = Month, y = 0 , yend = 10), #marking births
arrow = arrow(length = unit(0.25, "cm"), ends = "first"), size = 1)+
geom_text(data = births, aes(x = Month, y = 20, label = paste("Observed \n births")))+
geom_vline(aes(xintercept = "June"), size = 1)+#start of the monsoon ban
geom_vline(aes(xintercept = "August"), size = 1)+#end of the monsoon ban
geom_text(aes(x = "July", y = 80, label = "Monsoon ban"))+
scale_x_discrete(guide = guide_axis(n.dodge = 2))+
labs(y = "Snout to vent length (cm)")Proportion of neonates increases in March to May. Birth and neonates were observed around the same time.
Proportion of neonates in each month of 2018
snakes%>%
filter(Year == "2018", # data not sufficient for other years
Snout.to.Vent..cm. > 20)%>% # removing erroneous data
mutate(Month = factor(Month, levels = month.name))%>%
group_by(Species, Month)%>%
summarise(prop.neonate = sum(Snout.to.Vent..cm. < 40)/n())%>%
spread(Month, prop.neonate)| Species | January | February | March | April | May | October | November | December |
|---|---|---|---|---|---|---|---|---|
| Hydrophis curtus | 0 | 0 | 0.50 | NA | NA | 0 | 0.0769231 | 0 |
| Hydrophis schistosus | 0 | 0 | 0.02 | 0.05 | 0.0555556 | 0 | 0.0000000 | 0 |
Sex ratios
Proportion of females encountered by species
snakes%>%
group_by(Species)%>%
filter(Sex == "Male" | Sex == "Female")%>%
summarise(females = sum(Sex == "Female"),
N = n(),
P = females/N)| Species | females | N | P |
|---|---|---|---|
| Hydrophis curtus | 102 | 212 | 0.4811321 |
| Hydrophis schistosus | 304 | 618 | 0.4919094 |
Testing equality in sex ratios
# is proportion of females different from 0.5?
snakes%>%
filter(Sex == "Male" | Sex == "Female")%>%
group_by(Species)%>%
summarise(females = sum(Sex == "Female"),
N = n())%>%
group_by(Species)%>%
nest()%>%
mutate(test = map(data, ~prop.test(.$females, .$N, p = 0.5)), # Z - test for proportion
sumr = map(test, broom::tidy))%>%
select(sumr)%>%
unnest()%>%
select(-c(method, alternative))| Species | estimate | statistic | p.value | parameter | conf.low | conf.high |
|---|---|---|---|---|---|---|
| Hydrophis curtus | 0.4811321 | 0.2311321 | 0.6306857 | 1 | 0.4125066 | 0.5504525 |
| Hydrophis schistosus | 0.4919094 | 0.1310680 | 0.7173273 | 1 | 0.4518628 | 0.5320580 |
Sex ratio is not different from 0.5; p = 0.71.
Testing shifts in sex ratios over years
snakes%>%
filter(Sex == "Male" | Sex == "Female")%>%
group_by( Species, Year)%>%
summarise(N = n(),
females = sum(Sex == "Female"),
prop.female = females/N)%>%
ggplot(aes(Year, prop.female, fill = Species))+
geom_col(width = 0.5, col = "black", position = position_dodge())Logistic model to test for change in sex ratios of species across years:
snakes%>%
select(Species, Year, Sex)%>%
filter(Sex == "Male" | Sex == "Female")%>% # removing empty sex
mutate(Sex = factor(Sex, levels = c("Male", "Female")))%>%
group_by(Species)%>%
nest()%>%
mutate(mod = map(data, ~glm(Sex ~ Year, family = binomial, data = .)),
sumry = map(mod, ~tidy(.)),
r2 = map_dbl(mod, ~mcfadden(.)))%>%
select(sumry, r2)%>%
unnest(sumry)| Species | term | estimate | std.error | statistic | p.value | r2 |
|---|---|---|---|---|---|---|
| Hydrophis schistosus | (Intercept) | 0.5007753 | 0.2833779 | 1.7671643 | 0.0772007 | 0.0055997 |
| Hydrophis schistosus | Year2017 | -0.6120009 | 0.4379020 | -1.3975752 | 0.1622407 | 0.0055997 |
| Hydrophis schistosus | Year2018 | -0.6052544 | 0.2977861 | -2.0325137 | 0.0421017 | 0.0055997 |
| Hydrophis schistosus | Year2019 | -0.3404326 | 0.4010216 | -0.8489135 | 0.3959294 | 0.0055997 |
| Hydrophis curtus | (Intercept) | 1.0185696 | 0.3235748 | 3.1478639 | 0.0016447 | 0.0593767 |
| Hydrophis curtus | Year2017 | -1.3752445 | 0.5895404 | -2.3327403 | 0.0196618 | 0.0593767 |
| Hydrophis curtus | Year2018 | -1.4885732 | 0.3812161 | -3.9048020 | 0.0000943 | 0.0593767 |
| Hydrophis curtus | Year2019 | -1.2096248 | 0.4481189 | -2.6993392 | 0.0069477 | 0.0593767 |
Plotting predicted values:
snakes%>%
select(Species, Year, Sex)%>%
filter(Sex == "Male" | Sex == "Female")%>% # removing empty sex
mutate(Sex = factor(Sex, levels = c("Male", "Female")))%>%
group_by(Species)%>%
nest()%>%
mutate(mod = map(data, ~glm(Sex ~ Year, family = "binomial", data = .)),
fit = map(mod, ~predict(., years, se.fit = T, type = "response")),
fit = map(fit, ~data.frame(.)))%>%
select(fit)%>%
unnest(fit)%>%
mutate(Year = factor(c(2016:2019)))%>%
ggplot(aes(Year, fit))+
geom_pointrange(aes(ymin = fit - se.fit, ymax = fit + se.fit))+
labs(y = "P(Female|Year)")+
facet_wrap(~Species)+
theme(strip.text = element_text(face = "italic"))The proportion of female H. curtus changed significantly across years. H. schistosus sex ratios remained constant across sampled years.
Mortality in bycatch
Overall mortality rate
snakes%>%
group_by(Species)%>%
count(Condition.at.encounter..D.A.)%>%
# Calculating mortality rate
mutate(N = sum(n),
prop.dead = n/N)%>%
filter(Condition.at.encounter..D.A. == "Dead")%>%
select(-Condition.at.encounter..D.A.)| Species | n | N | prop.dead |
|---|---|---|---|
| Hydrophis curtus | 120 | 282 | 0.4255319 |
| Hydrophis schistosus | 174 | 967 | 0.1799380 |
Testing difference in mortality across species
Logistic model:
snakes%>%
filter(Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
select(Species, Condition.at.encounter..D.A.)%>%
mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")))%>%
nest()%>%
mutate(test = map(data, ~glm(Condition ~ Species, family = "binomial", data = .)),
sumr = map(test, ~broom::tidy(.)),
r2 = map(test, ~mcfadden(.)))%>%
select(sumr, r2)%>%
unnest()| term | estimate | std.error | statistic | p.value | r2 |
|---|---|---|---|---|---|
| (Intercept) | -0.2814125 | 0.1209241 | -2.327182 | 0.0199556 | 0.0502489 |
| SpeciesHydrophis schistosus | -1.2315652 | 0.1470904 | -8.372845 | 0.0000000 | 0.0502489 |
Fitted probability of mortality:
spp = data.frame(Species = c("Hydrophis schistosus", "Hydrophis curtus"))
snakes%>%
filter(Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
select(Species, Condition.at.encounter..D.A.)%>%
mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")))%>%
nest()%>%
mutate(mod = map(data, ~glm(Condition ~ Species, family = "binomial", data = .)),
fit = map(mod, ~predict(., spp, se.fit = T, type = "response")),
fit = map(fit, ~data.frame(.)))%>%
select(fit)%>%
unnest(fit)%>%
mutate(Species = factor(c("Hydrophis schsitosus", "Hydrophis curtus")))%>%
ggplot(aes(Species, fit))+
geom_pointrange(aes(ymin = fit - se.fit, ymax = fit + se.fit))+
labs(y = "P(Dead|Species)")+
theme(axis.text.x = element_text(face = "italic"))H. curtus had as significatly higher mortality rate in bycatch than H. schistosus
Mortality rate by age class
Are juveniles more vulnerable to bycatch mortality than adults?
snakes%>%
filter(!is.na(Class))%>%
group_by(Species, Class)%>%
count(Condition.at.encounter..D.A.)%>%
#calculating mortality rate
mutate(N = sum(n),
prop.dead = n/N)%>%
filter(Condition.at.encounter..D.A. == "Dead")%>%
select(-Condition.at.encounter..D.A.)| Species | Class | n | N | prop.dead |
|---|---|---|---|---|
| Hydrophis curtus | Adult | 40 | 62 | 0.6451613 |
| Hydrophis curtus | Juvenile | 41 | 125 | 0.3280000 |
| Hydrophis schistosus | Adult | 138 | 648 | 0.2129630 |
| Hydrophis schistosus | Juvenile | 5 | 105 | 0.0476190 |
| Hydrophis schistosus | Neonate | 3 | 8 | 0.3750000 |
Testing difference in mortality rates across age classes within species
Logistic model:
snakes%>%
filter(!is.na(Class),
Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
select(Species, Class, Condition.at.encounter..D.A.)%>%
mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")))%>%
group_by(Species)%>%
nest()%>%
mutate(test = map(data, ~glm(Condition ~ Class, family = "binomial", data = .)),
sumr = map(test, ~broom::tidy(.)),
r2 = map(test, ~mcfadden(.)))%>%
select(sumr, r2)%>%
unnest()| Species | term | estimate | std.error | statistic | p.value | r2 |
|---|---|---|---|---|---|---|
| Hydrophis schistosus | (Intercept) | -1.3012573 | 0.0960144 | -13.5527320 | 0.0000000 | 0.0300658 |
| Hydrophis schistosus | ClassJuvenile | -1.6944749 | 0.4681739 | -3.6193277 | 0.0002954 | 0.0300658 |
| Hydrophis schistosus | ClassNeonate | 0.7904317 | 0.7365814 | 1.0731085 | 0.2832225 | 0.0300658 |
| Hydrophis curtus | (Intercept) | 0.6931472 | 0.2738613 | 2.5310156 | 0.0113733 | 0.0823342 |
| Hydrophis curtus | ClassJuvenile | -1.3984157 | 0.3338240 | -4.1890813 | 0.0000280 | 0.0823342 |
| Hydrophis curtus | ClassNeonate | -16.2592154 | 1029.1215009 | -0.0157991 | 0.9873946 | 0.0823342 |
Fitted probability of mortality across age classes:
dev = data.frame(Class = c("Adult", "Juvenile", "Neonate"))
snakes%>%
filter(!is.na(Class),
Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
select(Species, Class, Condition.at.encounter..D.A.)%>%
mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")))%>%
group_by(Species)%>%
nest()%>%
mutate(mod = map(data, ~glm(Condition ~ Class, family = "binomial", data = .)),
fit = map(mod, ~predict(., dev, se.fit = T, type = "response")),
fit = map(fit, ~data.frame(.)))%>%
select(fit)%>%
unnest(fit)%>%
mutate(Class = factor(c("Adult", "Juvenile", "Neonate")))%>%
ggplot(aes(Class, fit))+
geom_pointrange(aes(ymin = fit - se.fit, ymax = fit + se.fit))+
labs(y = "P(Dead)")+
facet_wrap(~Species)+
theme(strip.text = element_text(face = "italic")) H. curtus adults had the highest mortality rate overall. H. schistosus juveniles had the lowest.
Mortality rate by sex
snakes%>%
filter(Sex != "")%>%
group_by(Species, Sex)%>%
count(Condition.at.encounter..D.A.)%>%
# Calculating mortality rate
mutate(N = sum(n),
prop.dead = n/N)%>%
filter(Condition.at.encounter..D.A. == "Dead")%>%
select(-Condition.at.encounter..D.A.)| Species | Sex | n | N | prop.dead |
|---|---|---|---|---|
| Hydrophis curtus | Female | 36 | 102 | 0.3529412 |
| Hydrophis curtus | Male | 52 | 110 | 0.4727273 |
| Hydrophis schistosus | Female | 62 | 304 | 0.2039474 |
| Hydrophis schistosus | Male | 72 | 314 | 0.2292994 |
Testing difference in mortality rates across sexes withing species
Logistic model:
snakes%>%
filter(Sex %in% c("Male", "Female"),
Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
select(Species, Sex, Condition.at.encounter..D.A.)%>%
mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")))%>%
group_by(Species)%>%
nest()%>%
mutate(test = map(data, ~glm(Condition ~ Sex, family = "binomial", data = .)),
sumr = map(test, ~broom::tidy(.)),
r2 = map(test, ~mcfadden(.)))%>%
select(sumr, r2)%>%
unnest()| Species | term | estimate | std.error | statistic | p.value | r2 |
|---|---|---|---|---|---|---|
| Hydrophis schistosus | (Intercept) | -1.3535045 | 0.1424629 | -9.5007493 | 0.0000000 | 0.0008565 |
| Hydrophis schistosus | SexMale | 0.1453737 | 0.1957905 | 0.7424963 | 0.4577867 | 0.0008565 |
| Hydrophis curtus | (Intercept) | -0.6061358 | 0.2071938 | -2.9254529 | 0.0034396 | 0.0133815 |
| Hydrophis curtus | SexMale | 0.5500463 | 0.2834464 | 1.9405655 | 0.0523110 | 0.0133815 |
Fitted probability of mortality across sexes:
sex = data.frame(Sex = c("Male", "Female"))
snakes%>%
filter(Sex %in% c("Male", "Female"),
Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
select(Species, Sex, Condition.at.encounter..D.A.)%>%
mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")))%>%
group_by(Species)%>%
nest()%>%
mutate(mod = map(data, ~glm(Condition ~ Sex, family = "binomial", data = .)),
fit = map(mod, ~predict(., sex, se.fit = T, type = "response")),
fit = map(fit, ~data.frame(.)))%>%
select(fit)%>%
unnest(fit)%>%
mutate(Sex = c("Male", "Female"))%>%
ggplot(aes(Sex, fit))+
geom_pointrange(aes(ymin = fit - se.fit, ymax = fit + se.fit))+
labs(y = "P(Dead)")+
facet_wrap(~Species)+
theme(strip.text = element_text(face = "italic")) There was no difference in mortality rates across sexes in either species.
Mortality and female reproductive status
Are gravid females more susceptible to bycatch mortality than the rest of our sample?
HS%>%
filter(Sex == "Female",
Class == "Adult")%>%
group_by(Gravid)%>%
count(Condition.at.encounter..D.A.)%>%
# Calculating mortality rate
mutate(N = sum(n),
prop.dead = n/N)%>%
filter(Condition.at.encounter..D.A. == "Dead")%>%
select(-Condition.at.encounter..D.A.)| Gravid | n | N | prop.dead |
|---|---|---|---|
| 41 | 151 | 0.2715232 | |
| Yes | 16 | 85 | 0.1882353 |
Testing thing difference in mortality rate of H. schistosus females by reproductive state
Logistic model:
HS%>%
filter(Sex == "Female",
Class == "Adult",
Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
select(Gravid, Condition.at.encounter..D.A.)%>%
mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")),
Gravid = ifelse(Gravid == "", "No", Gravid))%>%
nest()%>%
mutate(test = map(data, ~glm(Condition ~ Gravid, family = "binomial", data = .)),
sumr = map(test, ~broom::tidy(.)),
r2 = map(test, ~mcfadden(.)))%>%
select(sumr, r2)%>%
unnest()| term | estimate | std.error | statistic | p.value | r2 |
|---|---|---|---|---|---|
| (Intercept) | -0.9685592 | 0.1834379 | -5.280039 | 0.0000001 | 0.0087809 |
| GravidYes | -0.4929586 | 0.3326292 | -1.482006 | 0.1383386 | 0.0087809 |
Fitted probability of mortality by reproductive state:
gravid = data.frame(Gravid = c("Yes", "No"))
HS%>%
filter(Sex == "Female",
Class == "Adult",
Condition.at.encounter..D.A. %in% c("Dead", "Alive"))%>% # removing missing data
select(Gravid, Condition.at.encounter..D.A.)%>%
mutate(Condition = factor(Condition.at.encounter..D.A., levels = c("Alive", "Dead")),
Gravid = ifelse(Gravid == "", "No", Gravid))%>%
nest()%>%
mutate(mod = map(data, ~glm(Condition ~ Gravid, family = "binomial", data = .)),
fit = map(mod, ~predict(., gravid, se.fit = T, type = "response")),
fit = map(fit, ~data.frame(.)))%>%
select(fit)%>%
unnest(fit)%>%
mutate(State = c("Gravid", "Not Gravid"))%>%
ggplot(aes(State, fit))+
geom_pointrange(aes(ymin = fit - se.fit, ymax = fit + se.fit))+
labs(y = "P(Dead)")+
theme(strip.text = element_text(face = "italic"))Observed reproductive cycle of H. schistosus
Proportion of gravid females in the sample
# percentage of gravid females in sample
HS%>%
count(Gravid)%>%
mutate(N = sum(n))%>%
mutate(prop.gravid = n/N)%>%
filter(Gravid == "Yes")%>%
select(-c(Gravid))| n | N | prop.gravid |
|---|---|---|
| 88 | 967 | 0.0910031 |
Number of gravid females encountered per year
# checking the number of gravid females per year
HS%>%
group_by(Year)%>%
filter(Gravid == "Yes")%>%
count(Gravid)%>%
spread(Gravid, n)%>%
rename(`n(Gravid Females)` = Yes)| Year | n(Gravid Females) |
|---|---|
| 2016 | 1 |
| 2018 | 75 |
| 2019 | 12 |
Monthly variation in proportion of gravid females
Proper sampling was only conducted in 2018/19 and hence only this period is used for analysis of reproductive cycles.
# calculating the proportion of gravid females per month
gravid <- HS%>%
mutate(Month = factor(Month, levels = month.name))%>%
filter(Year == "2018" | Year == "2019")%>% # only for 2018/19
group_by(Month)%>%
summarise(N = n(),
gravid = sum(Gravid == "Yes"),
prop.gravid = gravid/N)
gravid| Month | N | gravid | prop.gravid |
|---|---|---|---|
| January | 102 | 12 | 0.1176471 |
| February | 134 | 34 | 0.2537313 |
| March | 129 | 24 | 0.1860465 |
| April | 54 | 7 | 0.1296296 |
| May | 36 | 1 | 0.0277778 |
| October | 16 | 0 | 0.0000000 |
| November | 69 | 3 | 0.0434783 |
| December | 34 | 6 | 0.1764706 |
Describing change in the relative proportions of gravid females across months and years of sampling.
# plotting prop gravid per month
gravid%>%
mutate(Month = factor(Month, levels = month.name))%>%
complete(Month, fill = list(prop.gravid = 0))%>%
ggplot(aes(Month, prop.gravid))+
geom_point(size = 3)+
geom_path(aes(group = 1), size = 1, linetype = "dashed")+
geom_segment(aes(x = "April", xend = "April", y = 0 , yend = 0.02), #marking births
arrow = arrow(length = unit(0.25, "cm"), ends = "first"), size = 1)+
geom_text(aes(x = "April", y = 0.04, label = paste("Observed \n births")))+
geom_vline(aes(xintercept = "June"), size = 1)+#start of the monsoon ban
geom_vline(aes(xintercept = "August"), size = 1)+#end of the monsoon ban
geom_text(aes(x = "July", y = 0.20, label = "Monsoon ban"))+
scale_x_discrete(guide = guide_axis(n.dodge = 2))+
labs(y = "Proportion of gravid females")Pregnancy from Novemeber to May. Peak in Feb. Observed two live births in April.
Development of embryos and eggs
Sample size
embryos%>%
summarise(n.mothers = length(unique(Field.Code)),
n.embryos = length(unique(Embryo.Code)))| n.mothers | n.embryos |
|---|---|
| 29 | 235 |
Summary statistics of egg measurements
embryos%>%
select(Egg.Length..mm.., Egg.Width..mm.., Egg.Weigth..g.., Snout.to.Vent..cm., Embryo.Weight..g.)%>%
skimr::skim()%>%
skimr::yank("numeric")%>%
select(c(skim_variable, mean, sd))Variable type: numeric
| skim_variable | mean | sd |
|---|---|---|
| Egg.Length..mm.. | 39.71 | 11.72 |
| Egg.Width..mm.. | 32.23 | 9.01 |
| Egg.Weigth..g.. | 14.30 | 5.45 |
| Snout.to.Vent..cm. | 17.72 | 5.24 |
| Embryo.Weight..g. | 4.83 | 6.28 |
Sex ratios within clutches
embryos%>%
group_by(Field.Code)%>%
filter(Sex != "")%>%
summarise(prop.female = sum(Sex == "Female")/n())%>%
skimr::skim(prop.female)%>%
skimr::yank("numeric")%>%
select(c(skim_variable, mean, sd))Variable type: numeric
| skim_variable | mean | sd |
|---|---|---|
| prop.female | 0.53 | 0.31 |
Number of feamle embryos
embryos%>%
filter(Sex != "")%>%
summarise(females = sum(Sex == "Female"),
N = n())| females | N |
|---|---|
| 86 | 166 |
Testing clutch sex ratio
broom::tidy(prop.test(86, 166, p = 0.5))%>% # Z - test
select(estimate:conf.high)| estimate | statistic | p.value | parameter | conf.low | conf.high |
|---|---|---|---|---|---|
| 0.5180723 | 0.1506024 | 0.6979603 | 1 | 0.4395567 | 0.5957383 |
The sex ratio in clutches is not significantly different from 0.5.
Infertility rate
Percentage of eggs with no embryos present:
embryos%>%
count(Embryo)%>%
mutate(N = sum(n), rate = n*100/N)%>%
select(-N)| Embryo | n | rate |
|---|---|---|
| Absent | 5 | 2.12766 |
| Present | 230 | 97.87234 |
Matrotrophy
Do female H. schistosus expend energy in the development of embryos? or Are the eggs formed and yolk only provides nourishment?
Change in yolk weight with embryo weight
embryos%>%
ggplot(aes(Embryo.Weight..g., Yolk.Weight..g.))+
geom_point(size = 3)+
geom_smooth(method = "lm", linetype = "dashed", size = 1)+
labs(x = "Embryo weight (g)", y = "Yolk weight (g)")Linear model to test change in yolk weight with embryo development:
embryos%>%
select(Yolk.Weight..g., Embryo.Weight..g.)%>%
nest()%>%
# Yolk weight is normally distributed
mutate(mod = map(data, ~lm(Yolk.Weight..g. ~ Embryo.Weight..g., data = .)),
sumr = map(mod, broom::tidy),
stat = map(mod, broom::glance))%>%
select(sumr, stat)%>%
unnest()%>%
select(term:r.squared)| term | estimate | std.error | statistic | p.value | r.squared |
|---|---|---|---|---|---|
| (Intercept) | 13.0905565 | 0.2965229 | 44.14687 | 0 | 0.4071222 |
| Embryo.Weight..g. | -0.4526283 | 0.0450509 | -10.04705 | 0 | 0.4071222 |
Yolk weight decreases by 0.45 g for every 1g increase in embryo weight.
Change in egg mass with embryo development
embryos%>%
ggplot(aes(Embryo.Weight..g., Egg.Weigth..g..))+
geom_point(aes(col = Yolk.Weight..g.), size = 3)+
geom_smooth(method = "lm", linetype = "dashed", size = 1)+
scale_color_viridis(name = "Yolk Weight (g)")+
labs(x = "Embryo Weight (g)", y = "Total Egg Weight (g)")Linear model to test relation ship between total egg mass and embryo development:
embryos%>%
select(Egg.Weigth..g.., Embryo.Weight..g.)%>%
nest()%>%
# Total egg mass is normally distributed
mutate(mod = map(data, ~lm(Egg.Weigth..g.. ~ Embryo.Weight..g., data = .)),
sumr = map(mod, broom::tidy),
stat = map(mod, broom::glance))%>%
select(sumr, stat)%>%
unnest()%>%
select(term:r.squared)| term | estimate | std.error | statistic | p.value | r.squared |
|---|---|---|---|---|---|
| (Intercept) | 13.5441434 | 0.4114909 | 32.91481 | 0 | 0.2188866 |
| Embryo.Weight..g. | 0.3676922 | 0.0557912 | 6.59050 | 0 | 0.2188866 |
Total egg mass increases by 0.36 g for every 1g increase in embryo weight.
Note: Embryo weight is log-normally distributed.
Along with decrease of fat bodies of gravid females (Figure A1) as embryos develop indicates matrotrophic nutrition.
Reproductive effort
Does the amount of enery expended by female H. schistosus reduce with age?
# Cleaning reproductive effort data
re <- embryos%>%
left_join(snakes, by = c("Date", "Field.Code", "Species"))%>%
select(Date, Field.Code, Embryo.Code, Species, Egg.Length..mm..:Tail.Length..cm..x,
Egg.Weigth..g..:Sex.x, Snout.to.Vent..cm..y:Tail.Length..cm..y, Weight..g.,
-Body.Length..cm..x, - Body.Length..cm..y)%>%
filter(Species == "Hydrophis schistosus")%>%
rename(Embryo.SVL = Snout.to.Vent..cm..x,
Embryo.TL = Tail.Length..cm..x,
Embryo.Sex = Sex.x,
Female.SVL = Snout.to.Vent..cm..y,
Female.TL = Tail.Length..cm..y,
Mother.Wt = Weight..g.
)%>%
# Removing erroneous data point
filter(Female.SVL > 50)Increase in clutch size with female age
re_clutch <- re%>%
select(Field.Code, Egg.Weigth..g.., Mother.Wt, Female.SVL)%>%
group_by(Field.Code)%>%
summarise(clutch.size = n(),
Clutch.wt = sum(Egg.Weigth..g..),
Mother.wt = last(Mother.Wt),
Female.SVL = last(Female.SVL))%>%
mutate(rcm = Clutch.wt/(Mother.wt - Clutch.wt))
re_clutch%>%
ggplot(aes(Female.SVL, clutch.size))+
geom_point(size =3)+
geom_smooth(method = 'glm', method.args = list(family = 'poisson'))+
labs(x = "Female SVL (cm)", y = "Clutch Size")Generalised linear model to test change in clutch size with female SVL:
Clutch size is poisson distributed with mean 8.1052632 and variance 10.3216374
re_clutch%>%
select(Female.SVL, clutch.size)%>%
nest()%>%
mutate(mod = map(data, ~glm(clutch.size ~ Female.SVL, data = ., family = "poisson")),
sumr = map(mod, broom::tidy),
stat = map(mod, broom::glance),
r.squared = map_dbl(stat, ~mcfadden(.)))%>%
select(sumr, stat,r.squared)%>%
unnest()%>%
select(c(term:p.value, r.squared))| term | estimate | std.error | statistic | p.value | r.squared |
|---|---|---|---|---|---|
| (Intercept) | -0.3764929 | 0.7174214 | -0.5247862 | 0.5997318 | 0.560094 |
| Female.SVL | 0.0256926 | 0.0073002 | 3.5194451 | 0.0004325 | 0.560094 |
The number of eggs borne by females increase by 1.0253151 for every 1 cm increase in female SVL.
Change in overall reproductive effort with age
Beta regression to determine relationship between relative clutch mass and female SVL:
library(betareg)
re_clutch%>%
select(Female.SVL, rcm)%>%
nest()%>%
mutate(mod = map(data, ~betareg(rcm ~ Female.SVL, data = .)),
sumr = map(mod, ~tidy(., conf.int = T)),
stat = map(mod, broom::glance))%>%
select(sumr, stat)%>%
unnest()%>%
select(component:pseudo.r.squared)| component | term | estimate | std.error | statistic | p.value | conf.low | conf.high | pseudo.r.squared |
|---|---|---|---|---|---|---|---|---|
| mean | (Intercept) | 0.6463332 | 1.6515322 | 0.3913537 | 0.6955358 | -2.5906105 | 3.8832769 | 0.0497975 |
| mean | Female.SVL | -0.0144713 | 0.0172664 | -0.8381195 | 0.4019636 | -0.0483129 | 0.0193702 | 0.0497975 |
| precision | (phi) | 12.8451280 | 4.7017084 | 2.7320129 | 0.0062949 | 3.6299489 | 22.0603071 | 0.0497975 |
re_mod <- betareg(rcm ~ Female.SVL, data = re_clutch)re_clutch%>%
mutate(fitted = predict(re_mod, re_clutch),
q.low = predict(re_mod, re_clutch, "quantile", at = c(0.05)),# lower bound
q.high = predict(re_mod, re_clutch, "quantile", at = c(0.95)))%>% # upper bound
ggplot(aes(Female.SVL, rcm))+
geom_point(aes(col = clutch.size), size = 3)+
geom_smooth(aes(y = fitted, ymin = q.low, ymax = q.high), stat = "identity", linetype = "dashed")+
scale_y_continuous(name = "Relative clutch mass")+
scale_x_continuous(limits = c(85, 115), name = "Female SVL (cm)")+
scale_color_viridis(name = "Clutch size")+
theme(legend.text = element_text(size = 10))The overall reproductive effort does not change with female age.
Change in reproductive effort per embryo with female age
Does the effort per embryo change with female age?
# calculating reproductive effort per embruo
re_embryo <- re%>%
select(Field.Code, Embryo.Code, Egg.Weigth..g.., Mother.Wt, Female.SVL, Embryo.Sex, Embryo.SVL)%>%
group_by(Field.Code)%>%
mutate(clutch.wt = sum(Egg.Weigth..g..),
Female.wt = Mother.Wt - clutch.wt)%>%
group_by(Field.Code, Embryo.Code)%>%
summarise(inv = Egg.Weigth..g../Female.wt,# investment per embryo
Female.SVL = last(Female.SVL),
Embryo.Sex = Embryo.Sex,
Embryo.SVL = Embryo.SVL)%>%
ungroup()Linear model to detemine relationship between female SVL and investment per embryo:
We have controlled for the effect of embryo development.
re_embryo%>%
select(Female.SVL, inv, Embryo.SVL)%>%
nest()%>%
mutate(mod = map(data, ~betareg(inv ~ Female.SVL + Embryo.SVL, data = .)),
sumr = map(mod, broom::tidy),
stat = map(mod, broom::glance))%>%
select(sumr,stat)%>%
unnest()%>%
select(component:p.value, pseudo.r.squared)| component | term | estimate | std.error | statistic | p.value | pseudo.r.squared |
|---|---|---|---|---|---|---|
| mean | (Intercept) | -2.2798667 | 0.4894683 | -4.657843 | 0.0000032 | 0.5627972 |
| mean | Female.SVL | -0.0161182 | 0.0050085 | -3.218169 | 0.0012901 | 0.5627972 |
| mean | Embryo.SVL | 0.0491682 | 0.0077622 | 6.334311 | 0.0000000 | 0.5627972 |
| precision | (phi) | 233.3499725 | 43.5079611 | 5.363386 | 0.0000001 | 0.5627972 |
Plotting relatonship between residuals and female SVL:
# using residuals to control for embryo development
emsvlinv <- betareg(inv ~ Embryo.SVL, data = re_embryo)
# Plotting
re_embryo%>%
modelr::add_residuals(emsvlinv)%>%
ggplot(aes(Female.SVL, resid))+
geom_point(aes(col = Embryo.SVL), size = 3)+
geom_smooth(method = "lm", linetype = "dashed", size = 1)+
scale_y_continuous(name = "Relative egg mass (residuals)")+
scale_x_continuous(limits = c(85, 115), name = "Female SVL (cm)")+
scale_color_viridis(name = "Embryo SVL (cm)")The relative egg mass (controlled for embryo development) reduces with female SVL.
Difference in reproductive effort for male and female embryos
Does female investment differ by sex of the embryo?
re_embryo%>%
filter(Embryo.Sex == "Male" | Embryo.Sex == "Female")%>%
modelr::add_residuals(emsvlinv)%>%
ggplot(aes(Female.SVL, resid, shape = Embryo.Sex))+
geom_point(aes(col = Embryo.SVL), size = 3)+
geom_smooth(method = "lm", linetype = "dashed", size = 1)+
scale_y_continuous(name = "Relative egg mass (residuals)")+
scale_x_continuous(limits = c(85, 115), name = "Female SVL (cm)")+
scale_color_viridis(name = "Embryo SVL (cm)")+
scale_shape_discrete(name = "Embryo Sex")Linear model to test for the difference in female investment by embryo sex:
We controlled for the effect of embryo development.
re_embryo%>%
filter(Embryo.Sex == "Male" | Embryo.Sex == "Female")%>%
select(Embryo.Sex, inv, Embryo.SVL)%>%
nest()%>%
mutate(mod = map(data, ~betareg(inv ~ Embryo.Sex + Embryo.SVL, data = .)),
sumr = map(mod, broom::tidy),
stat = map(mod, broom::glance))%>%
select(sumr, stat)%>%
unnest()%>%
select(component:p.value, pseudo.r.squared)| component | term | estimate | std.error | statistic | p.value | pseudo.r.squared |
|---|---|---|---|---|---|---|
| mean | (Intercept) | -3.7036341 | 0.1763377 | -21.003076 | 0.0000000 | 0.4781387 |
| mean | Embryo.SexMale | -0.1376065 | 0.0839749 | -1.638663 | 0.1012835 | 0.4781387 |
| mean | Embryo.SVL | 0.0481588 | 0.0083489 | 5.768291 | 0.0000000 | 0.4781387 |
| precision | (phi) | 205.9146533 | 38.4135707 | 5.360466 | 0.0000001 | 0.4781387 |
Female investment does not differ significantly with embryo sex.